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ABSTRACT 

We show that orbit-superposition dynamical models (Schwarzschild's method) provide reliable estimates of 
nuclear black hole masses and errors when constructed from adequate orbit libraries and kinematic data. We thus 
rebut two recent papers that argue that BH masses obtained from this method are unreliable. These papers claim to 
demonstrate that the range of allowable BH masses derived from a given dataset is artificially too narrow as a result 
of an inadequate number of orbits in the library used to construct dynamical models. This is an elementary error 
that is easily avoided. We describe a method to estimate the number and nature of orbits needed for the library. 
We provide an example that shows that this prescription is adequate, in the sense that the range of allowable BH 
masses is not artificially narrowed by use of too few orbits. This is illustrated by showing that the versus 
BH-mass curve does not change beyond a certain point as more orbits are added to the library. At that point, the 
phase-space coverage of the orbit library is good enough to estimate the BH mass, and the profile provides a 
reliable estimate of its errors. 

A second point raised by critics is that kinematic data are generally obtained with insufficient spatial resolution 
(compared to the BH radius of influence) to obtain a reliable mass. We make the distinction between unreliable 
determinations and imprecise ones. We show that there are several different properties of a kinematic dataset 
that can lead to imprecise BH determinations (insufficient resolution among them), but none of the attributes we 
have investigated leads to an unreliable determination. In short, the degree to which the BH radius of influence 
is resolved by spectroscopic observations is already reflected in the BH-mass error envelope, and is not a hidden 
source of error. The BH masses published by our group and the Leiden group are reliable. 

Subject headings: galaxies: nuclei — galaxies: statistics — galaxies: general 



L INTRODUCTION 

There are almost twenty detections of massive black holes 
(BHs) in galaxy centers that employ the technique of orbit su- 
perposition modeling (van der Marel etal. 1 998 ; Cretton & van 
den Bosch 1999; Cretton ef oL 1999; Gebhardt ef oL 2000a, 
2003; Cappeflari et al. 2002; Verolme et al. 2002b). The 
orbit superposition technique is based on a method originally 
invented by Schwarzschild (1979), who noted that the time- 
averaged orbits in a stationary potential can be summed (at fi- 
nite spatial resolution) to match the mass distribution that gives 
rise to the potential, thereby producing an equilibrium dynami- 
cal model. The Schwarzschild method was first used to analyze 
kinematic data for evidence of BHs in galaxy centers by Rich- 
stone & Tremaine (1985) for M87 and Dressier & Richstone 
(1988) for M31 and M32. Those early models were spheri- 
cally symmetric. In modern implementations of this method, 
large sets of orbits are run in a specified axisymmetric poten- 
tial (based on the starlight distribution and a central point mass) 
and then a non-negative linear superposition of orbits is found 



that best matches the kinematic and photometric observations. 
The usual goal of this process is the determination of only two 
quantities, the stellar mass-to-light ratio T (assumed to be in- 
dependent of position) and the BH mass M, . 

Valluri, Memtt & Emsellem (2004) (hereafter VME) and 
Cretton & Emsellem (2004) have challenged the reliability of 
BH mass determinations obtained via orbit superposition mod- 
eling. Cretton & Emsellem also argue that the uncertainties 
noted by VME can be dealt with via regularization (smoothing 
the distribution of orbit weights). 

VME raise a long list of problems and features of the orbit- 
superposition analysis technique. They focus on two main is- 
sues: orbit insufficiency and data insufficiency. First, they con- 
tend that the BH mass determination is artificially narrowed 
through the use of too small an orbit library and imply, but do 
not show, that we have made this mistake. In addition they sug- 
gest that the absence of a flat-bottomed x^ profile as a function 
of BH mass indicates that this problem has occurred. Second, 
they suggest that the BH mass M, can only be determined if 
the radius of influence r,- = GM,/a^ where the BH dominates 
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the stellar dynamics is well-resolved spectroscopically (here a 
is the line-of-sight velocity dispersion). In each of their main 
points they have fastened on one facet of rather complex prob- 
lems. 

We show in §2 that the use of an orbit library which is too 
small is an easily avoided elementary trap. In §3 we refer the 
reader to other work in which we show that we have adequate 
resolution to determine BH masses reliably, and argue that the 
presence of extensive radial kinematic coverage mitigates the 
uncertainties associated with a "not well-resolved" radius of in- 
fluence. We also show that the "flat-bottomed" behavior of 
is a feature whose presence depends on several features of the 
quality and extent of the kinematic data, not just the size of the 
orbit library. 

The principal result is that a profile constructed using our 
method does indicate the range of acceptable BH masses in a 
manner easily judged by the readers of our papers. 

2. ENOUGH IS enough: ORBIT SUFFICIENCY 

There appear to be at least three independent axisymmet- 
ric orbit-superposition codes, one developed by various Leiden 
students starting with van der Marel and continuing through 
Verolme and Cappellari, another developed by Gebhardt and 
Richstone and used by our collaboration (the "nukers") since 
1996, and a third described by VME. All practitioners agree 
that the orbits can be sampled by examining launch points in 
the three-dimensional space spanned by the isolating integrals: 
the energy per unit mass E, the z-component of the angular mo- 
mentum and the third integral 73. We further agree that this 
space is finite — only bound orbits matter (E < 0), the range of 

at fixed energy lies between zero and the angular momen- 
tum of a circular orbit of that energy, and the distribution of the 
third integral can be sampled by launching the orbits with zero 
velocity from positions spread across the zero velocity surface 
at fixed E and The remaining question is how densely or 
sparsely spaced these orbital parameters need to be (and there- 
fore how many orbits are required). The answer depends in part 
on the spatial resolution sought in the model, and in part on the 
complexity of the orbit structure in phase space, and therefore 
on the mass distribution. 

Nonetheless, it is possible to estimate the number of orbits 
required to represent all possibilities. The eventual construc- 
tion of the model requires matching the (stellar) mass density 
in a set of bins spanning the physical space of the (axisym- 
metric) model This bin set is composed of one radial dimen- 
sion (of tir bins, which are logarithmically spaced except very 
near the center) and one angular dimension (of np bins, which 
are equally spaced in sin(/5) where [3 is the latitude). To fully 
sample the equatorial orbits, we follow orbits with apocenter 
and pericenter in all possible pairs of radial bins; this requires 
roughly orbits. We then assign E and /, from the apo- 
peri sampling of the equatorial orbits and then use the same 
E and /; for non-equatorial orbits. For each E and L pair, 
we drop a set of orbits from points on the zero-velocity curve 
equally spaced in sin(/3) and more finely spaced than the an- 
gular bins. The models described here have 20 radial bins 
and 5 angular bins, hence our nominal orbit library contains 
Norb = 2 X «2/2 X X = 2 X 20^/2 x 25 =10,000 orbits. The 
initial factor of 2 comes from including each orbit's retrograde 
mirror image in the library, and the trailing factor of 25 results 
from oversampling the 5 angular bins by a factor of fp = 5. 

Even if one does not trust this estimate, it is straightfor- 



ward to investigate whether the orbit library is large enough 
by changing the number of orbits in the library. If the results 
are independent of the number of orbits then presumably the 
library is large enough. 
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Fig. 1. — The x~ goodness of fit obtained by comparing model kinemat- 
ics to observed kinematic data for NGC 891, versus model black hole mass. 
The models were constructed with successively finer coverage of phase space 
(and therefore larger numbers of orbits in the library), as described in the text. 
Once the number of orbits (in this example) exceeds 10,000 the x' profile is 
independent of the size of the library, except for small vertical offsets which 
we have removed for plotting purposes. In all these experiments we used the 
kinematic dataset described in §3. 

We created an orbit library as described above to model the 
galaxy NGC 821 and illustrate the results in Figure 1. The ob- 
servations are described in Pinkney et al. (2003) and in §3 
and the modeling is described in Gebhardt et al. (2003) We 
chose NGC 821 to correct an error in the models in Gebhardt 
et al. (2003), which used an incorrect Hubble Space Telescope 
(HST) point-spread function that did not reflect the spatial bin- 
ning of the STIS CCD (this was the only one of the 12 galaxies 
in that paper in which this error was made). The revised mass 
is (8.5 ±3.5) x 1O^M0, at a distance of 24.1 megaparsecs, a 
factor of 2.3 higher than was reported previously (at the same 
distance). In order to conduct the experiment illustrated in Fig- 
ure 1 we oversampled the grid by a factor of 2 relative to the 
prescription above in each dimension and then dropped out or- 
bits at random to obtain the libraries illustrated. In Figure 1 we 
show the x^ profile as a function of BH mass (we have akeady 
marginalized over the galaxy mass-to-light ratio). 

For small orbit libraries. Figure 1 exhibits the ragged x^ em- 
phasized by VME. It also shows that for our nominal library 
with Norb =10,000 orbits, the x^ profile is smooth enough to 
infer BH mass estimates and uncertainties, and that that larger 
libraries do not alter the mass estimate. Once the number of 
orbits is large enough that the x^ profile does not change as it is 
increased further, the inferred BH mass does not depend on the 
number of orbits, and the phase-space coverage must be good 
enough. 
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3. MORE IS BETTER: DATA SUFFICIENCY 

A second major point made by VME is that estimates of BH 
mass are unreliable if the BH radius of influence r, is not "well- 
resolved" by the kinematic data. This is one facet of the rich 
question of what different kinds of kinematic data reveal about 
the gravitational field of the galaxy. The geometry of the kine- 
matic data can be characterized by: (1) the spatial resolution 
compared to r,-; (2) the radial extent of the data (how far out the 
data extend); (3) the angular coverage (how many position an- 
gles, or whether there is integral-field data); and (4) the sparse- 
ness of the data (in radius and angle). In addition to these geo- 
metrical characteristics, the quality of the data (signal-to-noise 
and systematic errors such as template mismatches) also deter- 
mines what one can measure. A complete investigation of all 
these issues is beyond the scope of this paper. We emphasize 
the distinction between precision (are the error bars small?) and 
reliability (are the error bars accurately estimated?). We shall 
argue that geometric limitations to the quality of the kinematic 
data of the four kinds listed above reduce the precision of es- 
timates of M,, but the range of acceptable masses can still be 
judged from the profile, and therefore the mass determina- 
tion is reliable. 

We parameterize the resolution of the kinematic data near 
the center in terms of K = r,/6'D, where r, is the BH radius of 
influence defined earlier, D is the distance, and 9 is the full- 
width-half-maximum telescope resolution. VME argue that the 
resolution K must be much greater than unity for accurate BH 
mass determinations, although they do not give a clear numeri- 
cal criterion. A number of the BH mass determinations have K 
only slightly larger than unity. In the example of NGC 821 dis- 
cussed below, d ~ 0"08 (the observations were made at 8500A 
with a ."1 slit), D = 2A. \pc and r,- - %pc, so K - 0.9. 

VME's argument that observations with this resolution can- 
not determine M, is made without regard to the signal-to-noise 
of the data available. Even in the limit H = 0, excellent S/N data 
can reveal the presence of a BH. Given a model with constant 
mass-to-light ratio T, a central mass M, and a perfect LOS VD 
measured with H ^ 1, the shape of the LOSVD at low velocities 
(comparable to velocity dispersion of the galaxy as a whole) de- 
termines T, and the shape and extent of the high-velocity wings 
of the LOSVD determines the mass of a central BH. 

We have investigated the effect of varying the resolution 
K on our determinations of M, by comparing our masses 
determined using both high-resolution HST data and low- 
resolution ground-based data to masses determined using the 
same method with ground-based data alone. The addition of 
HST data generally improves K by a factor of 5, from less than 
unity to greater than unity. The results of these experiments 
are shown in Figure 8 in Gebhardt et al. (2003) and Figure 4 
in Kormendy (2004). We find that improved resolution gener- 
ally improves the precision of the measured M, by narrowing 
the profile, but the improved best-fit M, always lies within 
the range given by the estimated errors from the low-resolution 
data. This experiment is evidence that the profiles do exactly 
what they are supposed to — they provide an estimate of M, and 
of its uncertainty. Lower resolution data yield less precise val- 
ues of M,, but not unreliable ones. 

We next investigate the influence of kinematic data at larger 



radii on the precision of the BH measurement (points 2 and 3 
above), using data from NGC 821 (see Figure 2). This example 
shows that the x^ profile broadens and flattens as data at large 
radii and along the minor axis are discarded. The changes in the 
X^ profile imply a decrease in the precision of the measurement 
of M, and indicate the increased errors. If only HST data is 
used then models with and without a BH are both acceptable, 
because of the degeneracy between M, and T when only data 
inside and near the radius of influence are used. Note that the 
one model with only major axis data (the dashed line) yields 
a flat-bottomed x^ profile and a very uncertain estimate of M, 
compared to the model with additional data along the minor 
axis. We conclude that inadequate spatial coverage of kinemat- 
ics far from the BH can lead to a poorly defined BH mass and 
a flat-bottomed just as readily as can poorly resolved data 
(K > 1) in the galaxy center). Again, inadequate spatial cov- 
erage leads to imprecise BH mass estimates, but not unreliable 
ones. 
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Fig. 2. — Goodness of fit versus model black-hole mass in for differ- 
ent sets of observational data from NGC 821. All models were constructed 
with the 10.000 orbit library described in §2. All of these datasets comprise 
continuous (i.e. there are no gaps in spatial coverage) long-slit spectra from 
space- and ground-based observations. The solid lines illustrate the effect (and 
desirability) of extending the najor-axis data to larger radii. The dotted profile 
shows the effect of removing minor-axis data. If the errors in M, are deter- 
mined by the criterion Ax~ = 1 then all of the profile minima are consistent 
withM. = 8.5 X lO^M©. 

Turning to point 4, sparse radial sampling of the kinematic 
data (which we never do) can also produce a flat bottomed 
X^ profile [as illustrated in Figure 3 in Richstone & Tremaine 
(1985)] by pushing up or down velocity moments at unobserved 
radii to maximize or minimize model parameters. 

Kinematic data with better geometry (both resolution and 
spatial coverage) leads to smaller uncertainties in the BH mass. 
However, in all of the examples discussed above, the shape of 
the x^ profile provides an accurate estimate of this uncertainty, 
no matter how poor the data geometry may be. In the absence of 

'~ The complete set of observations modeled in Figure 2 consists of STIS spectra and ground-based data obtained at MDM Observatory. The STIS data were observed 
through the Of 1 wide slit with the central pixel centered on the galaxy and the data aggregated in a set of "observations" binned in rectangles with outer edges at 



0"025, Of'075, 0"125, 0"225, 0"425, 0"625, 0"925. The ground observations were obtained through a 1"0 slit binned into rectangles with outer edges at If 
I'.'l, 3'.'6, S'.'S, 9'.'2, U'.'e, 23'.'2, 36f'2, Sl'.'O. The minor axis data only extend to 23 '.'2. 
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good spatial coverage, data that resolves the radius of influence, 
or sufficiently high signal to noise, the BH mass will simply not 
be well determined, but in our experiments the true mass stiU 
lies within our error bounds. 

4. DISCUSSION 

We have made a point of distinguishing between imprecise 
(the error bars are large) measurements of the BH mass M, and 
unreliable ones (the error bars are wrong). 

Figure 1 shows that increasing the number of orbits beyond 
our standard prescription does not broaden the range of BH 
masses allowed by our method. Contrary to the assertion of 
VME, we must be using enough orbits. 

The issue of data sufficiency is far more complex. In §3 we il- 
lustrate (or cite) several ways that data insufficiency can lead to 
imprecise M, estimates. We also argue that these estimates are 
imprecise, not unreliable. Probably the best evidence for this 
is the fact that our M, estimates obtained with genuine orbit- 
superposition models lie within their own errors of repeat stud- 
ies with greatly improved resolution (Gebhardt et al. 2003)^^. 

Flat-bottomed profiles — a special case of imprecision — can 
be produced by kinematic data with limited radial or angular 
coverage, by sparsely sampled data, or by insufficient resolu- 
tion. A fourth way to produce a flat bottomed profile is the 
use of noiseless datasets constructed from two-integral models 
(Cretton & Emsellem 2004). This procedure creates models 
that match the data perfectly. Since axisymmetric models have 
more freedom than two-integral models, a range of three inte- 
gral distribution functions wiU — ^by construction — provide an 
exact match to the data. Since is bounded by zero, a flat 
bottom in the plot is unavoidable. 

Cretton & Emsellem (2004) argue that regularization 
(smoothing the distribution of orbit weights) reduces the un- 
certainty in M,. Our models use a form of regularization (max- 
imum entropy) but only as a numerical technique to accelerate 
convergence: we iterate our models while steadily reducing the 
weight of the entropy to the best-fit criterion until there is no 
further change in the best-fit model. Verolme et al. (2002b) 



study models with and without regularization for M32 and find 
no difference in the best fit BH mass. 

This paper has examined the two most serious recent criti- 
cisms of BH mass estimates from orbit superposition. We con- 
clude that the issues raised by these criticisms have not led us to 
report erroneous BH masses. There are a number of other tests 
of our procedure that also indicate that our method is reliable: 

1. our method has been successfully blind-tested against 
(spherical) Jaffe models containing BHs; 

2. our method yields the same masses as the Leiden code 
when applied to the same data (for NGC 821 and M32); 

3. BHs weighed by this method display the same distri- 
bution of residuals from the mass-velocity dispersion 
relation as BHs weighed via gas dynamics (Richstone 
2004). 

The tests described above are posted on the Nuker Team Web- 
site 

(http://www.noao.edu/noao/staff/lauer/nukerhtml). 

The BH masses published by our group (Gebhardt et al. 
2000a; Bower ef a/. 2001; Gebhardt ef a/. 2003; Siopis ef a/. 
2004) are obtained using an adequate number of orbits. The 
plots of x^(^») we publish provide the reader with an accurate 
assessment of the precision of the results. So far as we can teU, 
BH masses published from the Leiden group (van der Marel 
et al. 1998; Cretton & van den Bosch 1999; Verolme et al. 
2002b; Cappellari et al. 2002) also should be reUable. 

We thank the director and staff of the Carnegie Institution 
Observatories for hospitahty during a meeting at which this pa- 
per was developed. Support for proposals 8591, 9106 and 9107 
was provided by a grant from the Space Telescope Science In- 
stitute, which is operated by the Association of Universities 
for Research in Astronomy, Inc, under NASA contract NAS 5- 
26555. This research was also supported by NASA Grant NAG 
5-8238. 
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^ We have not investigated filled kinematic maps in radius and angle, as might be obtained with integral-field devices like SAURON. Especially where these data 
imply triaxiality (thereby violating the basic assumption of our axisyimnetric models), there may be unexpected changes in M,. 



